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We consider heat transfer between two thermal reservoirs mediated by a quantum system using the gener- 
alized quantum Langevin equation. The thermal reservoirs are treated as ensembles of oscillators within the 
framework of the Drude-Ullersma model. General expressions for the heat current and thermal conductance are 
obtained for arbitrary coupling strength between the reservoirs and the mediator and for different temperature 
regimes. As an application of these results we discuss the origin of Fourier's law in a chain of large, but finite 
subsystems coupled to each other by the quantum mediators. We also address a question of anomalously large 
heat current between the STM tip and substrate found in a recent experiment. The question of minimum thermal 
conductivity is revisited in the framework of scaling theory as a potential application of the developed approach. 

PACS numbers: 05.70.Ln, 05.10.Gg, 65.80.-g 



I. INTRODUCTION 

The development of non-equilibrium thermodynamics 
mainly goes along two distinct directions. One direction is 
the study of the energy transfer through microscopic systems, 
such as nanotubes, molecules, or quantum dots Qyl. Beyond 
a purely academic interest in the problem, research suggests 
that nanoscale and molecular systems are good candidates 
for many technological advances, such as molecular wires, 
molecular diodes, rectifiers, and switches (HQ]. The other 
direction, with longer history, deals with thermalization of, 
and energy flow through finite, macroscopic subsystems. Ex- 
amples of such an approach are the Caldeira-Legett JUl and 
Nieuwenhuizen-Allahverdian j(| 0] models in which a ther- 
mal reservoir or macroscopic system is described as a large or 
infinite ensemble of harmonic modes. 

In microscopic systems, such as chains of multilevel sys- 
tems |Hl U, harmonic oscillators Il0l - [l2"ll . or spins lfl3ll the 
relaxation processes of individual elements and the processes 
of mutual equilibration between different elements are insep- 
arable from each other and all take place on a microscopic 
time scale. In contrast, the local equilibrium requirement for 
macroscopic systems implies that the equilibration processes 
proceed on two vastly different time scales. The local equi- 
librium is established on the microscopic time scale, while 
the equilibration between the macroscopic subsystems takes a 
much longer time. 

One of the most visible problems of non-equilibrium ther- 
modynamics is the microscopic derivation of Fourier's law, 
specifying that the heat flux j through both fluids and solids 
is given by j = -KVT(r), where the temperature T varies 
smoothly on microscopic scale and k is the thermal conduc- 
tivity. Despite the ubiquitous occurrence of this phenomenon, 
very few rigorous mathematical derivations of this law are 
known fTill . While for 3D generic models Fourier's law is 
expected to be true, this law may not be valid for ID and 2D 
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systems fl]. 

A recently developed approach to study heat transport at 
the microscopic level is based on either the classical or the 
quantum Langevin equation. The quantum Langevin equation 
was first considered in Ref. lfl5ll for a weakly damped har- 
monic oscillator. Later fltjll . it was used to formulate trans- 
port, collective motion, and the Brownian motion from a uni- 
fied, statistical-mechanical point of view. In Refs. JU-Q, fl7h - 
[l9ll . the Langevin equation was used for studying the thermal- 
ization of a particle coupled harmonically to a thermal reser- 
voir and other closely-related problems. This approach was 
generalized in Refs. Bill |20h22I1 in order to explore the non- 
equilibrium steady-state heat current and temperature profiles 
in chains of harmonic oscillators placed between two thermal 
baths. Recently, a new method for an exact solution to the 
Lindblad and Redfield master equations for open quadratic 
system of n fermions in terms of diagonalization of a An x An 
matrix has been developed 1I23II24I1 . The method has been ap- 
plied to Heisenberg XY spin 1/2 chain coupled to heat baths 
at its ends. Generally, this approach can be considered as an 
alternative to the quantum Langevin equation. 

In this paper we investigate the non-equilibrium steady- 
state heat transfer between two thermal reservoirs described as 
ensembles of harmonic modes mediated by a quantum system, 
which is also considered in the harmonic approximation. This 
is a Hamiltonian system with fixed total energy, but increas- 
ing entropy. Our approach is based on the quantum Langevin 
equation and uses the Drude-Ullersma model (DUM) for the 
bath mode spectrum. This is a generalization for the non- 
equilibrium situation of the approach employed in Ref. 0] to 
the study of statistical thermodynamics of a quantum Brown- 
ian particle coupled to a single thermal reservoir. 

The solution obtained within this model allows us to de- 
termine the heat conductance between two thermal baths at 
arbitrary strength of the coupling constant. The results pre- 
sented here are valid for an arbitrary temperature difference 
and arbitrary cut-off frequency, which plays the role of the 
Debye frequency. As is found, temperature dependence of the 
conductance may possess a plateau at intermediate tempera- 
ture range, similar to the "classical" plateau at high tempera- 
tures. Dependence of the thermal conductance on the coupling 



strength displays a maximum. We also show that the quan- 
tum thermal bath approach, in which a many-body problem is 
replaced by the one-body approximation where the effect of 
thermal baths is quantified by random forces, yields the same 
results for conductance as the rigorous (many-body) solution 
in the limit of large Debye frequency. On a more general note, 
this approximation can be successfully used for solving more 
complex problems without necessarily resorting on the rigor- 
ous solution that is based on the full-fledged Hamiltonian. 

The solution to the problem of heat transfer between two 
thermal reservoirs with a quantum particle as the mediator is 
applied to a chain system consisting of macroscopic subsys- 
tems coupled to each other by quantum particles (mediators). 
The microscopic time scale r describes the time it takes for the 
heat current facilitated by the mediators to come to a steady- 
state. Each subsystem has arbitrary large heat capacity and 
the equilibration time between them is much longer than r. 
In this case, Fourier's law follows naturally as the differential 
form of the energy conservation law. 

We use these results to explain a recent experiment in which 
the heat flow in vacuum between an STM tip and a substrate 
was found to be about ten orders of magnitude greater than 
that expected from the blackbody radiation theory ll25ll . Our 
suggestion is that the heat flow in this experiment was me- 
diated by a carbon monoxide molecule placed in the gap be- 
tween the tip and substrate. In addition, we briefly discuss 
the problem of minimum thermal conductivity attained when 
the coherence length of the phonons is minimal and of the or- 
der of the interatomic distance. Finally, we also briefly men- 
tion a possible application of the developed model to study 
the Josephson junction, which provides an important example 
of strong coupling between the quantum system and thermal 
baths. 

The paper is organized as follows. The model is introduced 
in Sec. II, where the generalized Langevin equation is derived 
and solved using the DUM. In Sec. Ill, expressions for the 
heat current between the thermal baths and heat conductance 
are derived and analyzed for different temperature regimes 
and different coupling strengths. In Sec. IV, we compare the 
solution obtained in the quantum thermal bath approach with 
the the rigorous solution. In Section V the specifics of strong 
coupling is described. Section VI is devoted to Fourier's law 
in a chain of macroscopic subsystems. Section VII discusses 
the application of our model to anomalously large heat flow 
between the STM tip and substrate. Sections VIII and IX 
discuss possible applications of the model to deal with the 
problem of minimum thermal conductivity and the Josephson 
junction, respectively. 
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FIG. 1: (Color online) Diagram representation of the total Hamil- 
tonian {TJ . The large circles correspond to the Hamiltonians of the 
thermal reservoirs, Eqs. (TJJ, the small circle is the Hamiltonian (fJJ of 
the mediator, and the dotted lines correspond to the coupling Hamil- 
tonians ([4} . 
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is the Hamiltonian of the quantum system (the mediator) de- 
scribed as a harmonic oscillator, 
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are the Hamiltonians of the first (v 
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describe interaction between the mediator and the baths. In 
Eq. (O, x and p are the coordinate and momentum operators 
and m and k are the mass and spring constant of the particle. 
In Eqs. (O and x vi and p vi are the coordinate and mo- 
mentum operators, whereas m v i and uj U i are the masses and 
frequencies of the oscillators for the ith mode that belong to 
the t4h bath. Finally, C V i are the coupling coefficients that 
describe the interaction between the quantum system and the 
baths. The last contributions to the right-hand side of are 
self-interaction terms, which guarantee that H^+V^ are pos- 
itively defined operators. Fig. Q]contains graphical represen- 
tation of the Hamiltonian (Q]). Here the large circles represent 
the Hamiltonians of the baths (01, the smaller central circle 
stands for the particle Hamiltonian @, and the dashed lines 
describe the interaction between the central particle and the 
baths ©. 

The Heisenberg equations for the coordinate and momen- 
tum operators corresponding to each bath and for our quantum 
system read 



II. LANGEVIN EQUATION 
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The total Hamiltonian of the system under consideration is 
similar to that in Refs. UHI^ 



Put ^^vi^iiiXvi ~t~ CviX : 
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Htot =H+H m + H B 2 + Vi + V 2 . 
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and 
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Considering Eqs. (O and |6) as a system of inhomogeneous 
equations with known inhomogeneity, its solution can be writ- 
ten as 

x vi (t) = x v i(Q) cos(w„ji) + Pvi ^ sm(uj vi t) + 



m vi uj vl 



C v 



m vi LO 
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ds sm[ij v i(t — s)]x(s) (9) 



and 



p„i{0) cos(u vi t) + C vi / dscos[u}vi(t— s)]x(s). (10) 
Jo 

Substituting (O into <[8J an d integrating by parts, one obtains 
the quantum Langevin equation: 



mx = -kx(t) + i](t)- / dsj(t-s)x(s)-j(t)x(0), (11) 
Jo 



where 



2^(0) cos(uuit) + Put ^ sm.{ujuit) 



(12) 



where k = 1,2,3, ...N v . In Eq. ( TBI , are the mode 
spacing constants, D v are the characteristic cutoff frequen- 
cies qualitatively similar to the Debye frequency, and j v are 
the coupling constants between a given reservoir and the me- 
diator 13]. Hereafter we assume for simplicity that D\ = 
D2 = D. In the final results we take the limit N v 00 
and A„ — > 0. 

Substituting Eq. dT4T > into Eq. ( fT3l ) and replacing summa- 
tion over i by integration, one arrives at the following expres- 
sion for 7(i): 



7(i) = "fDe 



-D\t\ 
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where 7 = 71+72- Using dT4b and (fT3T >. Eq. dTTb can be 
solved by Laplace transformation l29ll : 

x(t)=g(t)x{0) + -g(t)p(0) + - [ dsg(t- s) v (s). (16) 



Here 
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where g(z) is the Laplace transform L of g(t), L 1 is the 
inverse Laplace transform, and 



m L> + z to 
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After substituting ( fT8l into ( fTTI ). g(t) can be presented as 



is the noise that comes from the baths and 
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is the friction kernel, which takes into account the interaction 
of the quantum particle with both thermal reservoirs. 

A. Drude-Ullersma model 

At this point we have made no specific assumptions about 
the properties of the Hamiltonians that describe the thermal 
reservoirs. The microscopic structure of the thermal reser- 
voirs does not affect the nature of the energy exchange be- 
tween them. Therefore, one can choose a specific, physically 
meaningful model of the Hamiltonians Hbi and Hb2 without 
sacrificing the generality of the results. 

Here we employ the Drude-Ullersma model (DUM) |0,H3 
l28ll . The model assumes that in the absence of the interac- 
tion with the quantum system, each bath consists of uniformly 
spaced modes and introduces the following w-dependence for 
the coupling coefficients: 



uj vk = kA u , C vi 



1 2-y v m vi cjl i A v Dl 
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Here g n are defined by the last two relations and pi n are the 
roots of equation 



(/i-D)(M 2 +^)+7Au = 0, 



(21) 



where ljq = \J k/m. Statistical thermodynamics of a quan- 
tum particle coupled to a thermal bath was considered in J3l 
in the limit of large D when 

D » u> , 1/t p = 7/TO = 7, 1/t x = fc/7 = WqT p . (22) 

In this case, ( fT9l and (f2Tb give 

Hi,2 » ^— (iTr), [13 « D - l/r p » |/xi )2 | (23) 



and 
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where r = ^t p /t x . If t p /t x = (w t p ) 2 < 1/4, fj,i, 2 

are real and if t p << t x , the quantities r p and t x can be inter- 
preted as the characteristic relaxation times for the momentum 
and coordinate, respectively. Otherwise, when wot p > 1/2, 
/ii ; 2 = l/2T p =Fi/r , where t = (Wq — 7 2 /4) -1 / 2 determine 
the oscillation time while t p again determines the damping 
time. In what follows, however, we consider a more general 
case when D, ujq, and r" 1 may be comparable. In this case, 
the roots \i n and coefficients g n , where n = 1, 2, 3, are deter- 
mined as analytical solutions of ( fJTT i and from relation dT9b , 
respectively. 



III. HEAT CURRENT 

Using (0 and ||6), one can easily show that the rate of 
change of the energy of a given thermal reservoir is given by 



d / Pti , muiuiiXZi \ _ 



where 



2m,,, ; 



(26) 



is the work the quantum system performs on the ^th bath per 
unit of time (the power dispersed in the j/th bath) ll20ll . In 
the steady-state regime the power acquired by one reservoir is 
taken from the other, so that the steady-state heat current J t h 
can be presented as J t h = (Vi) = -(P2) or m the symmetrical 
form 



(27) 



Using solutions fox p U i(t) and x(t) fromEqs. ( TTOb and (fl6l l. 
and omitting here the transient processes that wash out over 
the short time r = max(T p , r x ) we get 



2m f— ' m„i 
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where 



£f = I dsg(t - S )(p m (0)v(s) + v(s)Pui(0)), (29) 
Jo 



$= dsg(t- s)^ (0)77(5) + v (s)x vi (0)), (30) 
Jo 



and 



jii = ( dr cosu vi (t - t)x(t) / dsg(t - s)<q(s) + 
Jo Jo 

/ dsg(t — s)rj(s) / dr cosw^i — t)x(t)). (31) 
Jo Jo 



The derivation of the steady-state expressions for j[ a \ j[ b \ 

(c) 

j\ and ultimately for Jth depends on how the contact be- 
tween the baths is established. A physically meaningful model 
should yield the same result regardless of how the coupling 
initially takes place. To verify this we have considered two 
options. In one case, the quantum system is attached simulta- 
neously to both baths at time t = 0. In the second case, the 
quantum particle is coupled first to the first bath, reaches ther- 
mal equilibrium with it, and at a later moment (which is again 
£ = 0) it is coupled to the second bath. 

We have established that the steady-state heat current is the 
same in either scenario. Below we give a derivation in the 
case of simultaneous coupling of the mediator to two thermal 
reservoirs at t = 0. We can assume that for t < the dynamic 
variables of the baths are determined by the usual expressions: 



Xvi (t) 



2m„,u l 



-iuj u it 



( 25 ) and 



p V i{t) = rn vi x vi {t). 



(32) 
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Here the creation and annihilation operators a„; and sat- 
isfy [a^i,a^i k ] = 5ik$w'- The operators' Gibbsian ensemble 
averages are determined by 



(ati^'k + aiy'ka^) = 5 ik 5 yv > coth ( 



where [i v = (/cbT^) 1 , which also result in 

(p«(0)jv,(0)) = m 2 ^ 2 ,^ (0)^,(0)) 
Oij<W coth I — — 



(34) 



and 



(p«(0)av,(0) + v 3 (0K,(0)) = 0. 



(35) 



(36) 



Using these equations, as well as Eq. 1121 . the ensemble av- 
erages (p vi (0)r)(t) + v{t)Pui(0)), (x ui (0)ri(t) + r)(t)x vi (0)) 
can be found and l(P v ) can be written as 
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where 
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Evaluating the integrals in d38l l using $1% and omit- 
ting exponentially decaying contributions, one can find 
that all time-varying terms, such as those proportional to 
sin(w,/,£) cos(w„ji), etc. cancel each other out. It means 
that the steady-state heat current is truly time-independent and 
does not contain any fluctuating contributions. Using Eq. ( [Pil l 
and substituting integration for summation we obtain 



(TV 



in 



h lv D 2 



duJUJ coth(/3 !y fio;/2) 
— - ■ "J (w* + D*)(u* + ft) 



(40) 



Similarly, Eq. (139b can be rewritten as follows 



(2) 



2 rt 



■^-(J o drS{t-r)x{r) J dsg(t-sHs) 
+ / dsg(t — s)rj(s) / dr cos S(t — t)x(t)}, 



where 



= E 



cos(w„t) = 2S(t) - De- Dt . (41) 



^ "f„D 2 m m 



Using ( [T2l , d35l l. and (1361 . (V u ) ^ can be obtained in a similar 
way, reaching its steady-state value when t —> oo, and the heat 
current is given by 



-/th 



HP 2 

2lTT r 



E 5n/X ™ 



dujui[ni(uj) — n 2 (u>)] 





(42) 



where n u (u>) = l/[cxp(fio;/3^) — l] are the phonon occupation 
number for the respective thermal reservoir. Here we assume 
for simplicity that 71 = 72- In a more general case when 
7i 7^ 72> ^th will be determined by the same expression d42l 
by substituting l/2r p — > 27172/(7771). It can be shown ex- 
plicitly that the same result for (V^) 1 - 2 ^ and, eventually, for 

(c) 

J t h can be obtained if one performs r-integration in j wi first 
and, finally, calculates the i-sum in (f39b . This provides an ad- 
ditional verification of the formula d4Tb and expression d42l . 
If |Ti - T 2 | < (Ti + T 2 )/2 = T, Eq. g2]> determines the 
heat conductance K: 



modes of the Hamiltonian "Hbi are shifted due to the interac- 
tion with the mediator to the values 



v\k = Uik - —4>{uk), 

7T 



(44) 



where <p(oj) is a certain known function of the parameters of 
the Hamiltonian. In the limit of small 71, 4>{uS) ~ 71. 

In the Appendix we have shown that after coupling of this 
thermalized combined system to the second thermal bath with 
a different temperature, the same steady-state heat current (1421 
will be established despite the small difference in the micro- 
scopic makeup of the two baths. This is important not only 
from the physical point of view that the steady-state energy 
current between two thermal reservoirs should not depend on 
the initial conditions, but also for the derivation of Fourier's 
law as shown below. 

It should be noted that the existence of a unique steady- 
state independent of the initial conditions cannot be taken for 
granted. There is substantial literature devoted to this prob- 
lem in classical and quantum systems HI l30l - l34ll . Ref. ll34ll 
considers the existence of the steady-state and thermal equi- 
libration in a system that represents a quantum wire coupled 
to two baths. It appears that a necessary, but not sufficient 
condition of uniqueness of the steady-state is absence of the 
bound state in the spectrum of the thermal bath. In the model 
we are using here the bound state would manifest itself as an 
imaginary root in Eq. (1211 . This would result in non-decaying 
oscillatory contribution to the steady-state energy current and 
its dependence on initial conditions. However, the absence of 
the bound states does not solely guarantee the uniqueness of 
the steady-state or thermal equilibration [01 [Hill. In this 
regard, our results demonstrate explicitly the existence and 
uniqueness of the steady-state in the considered model. 



A. Different temperature regimes 

Expressions for the heat current and heat conductance can 
be simplified in the limits of high and low temperatures when 
h\\i n \jk^T v <§C 1 and fi\[i n \l^'BT l , ^> 1, respectively. Here 
Mn in = 1,2, 3) are the roots of Eq. d2U . 

In the high temperature limit, Eqs. (1421 and (l43l l reduce to 
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AT^O AT 



rlk B D 2 



8ttt„ 



dww 2 cosech 2 ((3huj /2) 

(^ +W 2 )(M 2 +W 2) - 



(43) 



where AT = T 2 - T t and r h = h/k B T. 

In the second scenario of consecutive coupling of the me- 
diator to the thermal baths, one can consider initially a closed 
system that describes the equilibrium (Gibbsian) state of the 
first thermal bath plus the mediator. The corresponding set 
of eigenvalues and eigenmodes {vk,ek}, can be determined 
by diagonalization of the Hamiltonian H + Hbi + Vi. As 
was found in Ref. 0], the frequencies uiik of the unperturbed 



Jt h = -K(T 2 - Tl ); Kk-^Y,J^-. 



4t p ^ D + fj, r , 



The above sum can be written as 

3 



E 



{ D + fJ, r 



= -L[g(t)h=D = -Dg(D) (46) 



when using the well known properties of the Laplace trans- 
form and also relation J^n=i 5»Mri = [(~ l) fc ~ where 
k = 0, 1,2. Thus, K in ( f45l > can be written as 

fc R 2D 2 



K 



At p [2{D 2 +u 2 ) + ^D\ 



(47) 
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In the limit (1221 of large D we obtain 

k B fc B 7 



K ~t=*- (48) 

In the deep quantum regime (low temperatures), we find 

(49) 



^ ^kj{Tt-T%) ^g n 

Jth ^5 2^ 7^ 



30^ 3 r n 



n=l ^" 



Using again the Laplace transform and ( l20b . we find 



E 

n=l 



^g(^ = o) = j_ 

dz Wn 



(50) 



and finally have, in the quantum regime, 
J**JU^(T}-'I*) and KnfL^. (51) 



30/l 3 W 4 T 2 



15?i 3 w 4 t 2 ' 



The temperature dependence of Jth is the same as in the 
Stefan-Boltzmann law 



7T 2 fc 4 

JsB = A m^ {T ^ T ^ 



(52) 



where A is the area of a black body from which radiation 
emits. We will use this observation in Sec. IVIII for discus- 
sion of experimental results found in l25ll . 



IV. QUANTUM THERMAL BATHS 

Let us show now that the result for the heat current given 
by Eq. ( l42l in the limit D — > oo can be reproduced within the 
quantum thermal bath (QTB) approach. In the phenomeno- 
logical QTB model the many-body problem described by the 
Hamiltonian ([TJ is replaced by a one body problem in which 
friction is introduced "by hand", instead of being a logical 
consequence of energy redistribution between the practically 
infinite number of modes. The other effect of a thermal bath 
is modeled by a random force. The equation of motion for the 
mediator takes form of the Langevin equation 



mx + yx + mcjQX = F\ (t) + F 2 (t), 



(53) 



where the stochastic forces (colored noise) F\ and F 2 describe 
the effects of the two heat baths with temperatures T\ and 
T 2 , respectively. We take (F„(t)) = and (F v (t)F v .(t')) = 
S vu 'K v (t — t'), where v, v' = 1,2 and K„(t) is determined 
by its Fourier transform as 



*Mt> = L 



dujK v (u))e 



with 



7 



twcoth.(hu/2knT v ). 



(54) 



(55) 



The correlators of the Fourier transforms of the random forces 
are given by 

(F v (u))F v ,(oj')) = 2-k8{lo + u')8 vv ,K v {u). (56) 
The solution of Eq. d53l is given by 

1 F 1 (w) + F 2 (w) 



X(0J 



m oj 2 + i^lo — lJq 



, 7 = 7/m. (57) 



The energy conservation law for the mediator is given by 



= - 7 (± 2 ) + (F lX ) + (F 2 x), (58) 



dt \ 2 



which in the steady-state corresponds to j(x 2 ) = (Fix) + 
(F 2 x). 

We define heat current J$ as the energy transferred per 



'th 

unit of time from the first to the second bath: J, 



—dEi/dt, or 



(s) 
th 



= -\l(i 2 ) + (Fix) = i«Fi±) - (F 2 x)). (59) 

Here E\ is the internal energy of the first bath and we take 
that the energy dissipated by the mediator is equally split be- 
tween the two baths. In general, the energy dissipated by the 
mediator can be split between the thermal baths in arbitrary 
proportion, but then the correlators ( 1581 also will be propor- 
tional to the fraction of energy dissipated by the mediator in a 
given bath. 

Using Eqs. d57l i and d56l l. we obtain 



'th 



Atttyi 



duJU)[Kx(u) - K 2 (w)} 
u) 2 + i^w — uOq 



(60) 



Taking into account expression 031 ), the heat current ( f60b can 
be rewritten as 



'th 



27TT 2 



du>UJ 3 [ni(ui) — n 2 (ui)] 



(w 2 



a\2 



+ J 2 U! 2 



and the corresponding conductance is 



K (s) 



87TT 2 



dww 4 cosech 2 (fiftw /2) 
(uj 2 — cj 2 ) 2 + 7 2 w 2 



(61) 



(62) 



Thus, the heat current given by Eq. doTl ) and obtained from 
the phenomenological quantum thermal bath approach coin- 
cides in the limit D — > oo with that given by Eq. (l42l , which 
is the result of rigorous solution of a microscopic many-body 
model of the thermal reservoirs. This is a solid indication that 
the quantum thermal baths approach can be used to address 
more complicated problems of energy transfer via quantum 
mediators without resorting to a full-blown treatment based 
on a many-body Hamiltonian such as the one given by Eq. 
([TJ. It should be noted that the QTB model has been recently 
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successfully used by Dammak et al. 13511 for sampling quan- 
tum fluctuations within the framework of molecular dynam- 
ics (MD) simulations. Using the QTB model, the authors re- 
produced several experimental data at low temperatures in a 
regime where quantum statistical effects cannot be neglected. 
Our result here suggests that the MD approach can account 
for quantum statistical effects in non-equilibrium situations as 
well. 

Also we should mention that a powerful method for solv- 
ing problems involving open quadratic systems for fermions 
has been developed by Prosen iHl and Prosen and Zuncovic 
l24ll . It would be interesting to see how a problem of non- 
equilibrium bosonic systems, like the one considered here, 
can be reformulated in terms of this novel approach of "third 
quantization". 



V. WEAK AND STRONG COUPLING REGIMES 

In this section we return to the analysis of our results em- 
phasizing the effects of the weak and strong coupling on heat 
transfer as well as some interesting features in the behavior of 
the heat conductance. The main purpose of this section is to 
clarify the conditions which may allow us to assign a certain 
temperature to the mediator. This is possible when the media- 
tor is in a state of weak non-equilibrium. For a simple system 
with two degrees of freedom the weak non-equilibrium means 
that the virial theorem is approximately satisfied. This condi- 
tion depends on the coupling strength 7 and temperature. Let 
us consider the average potential and kinetic energy: 



&bTx = (kx 2 ) and 



(p 2 /m) 



(63) 



The imbalance between them has been proven to be useful in 
determining the statistics of a quantum particle coupled to a 
single heat bath 0]. Even when the total system - (heat bath + 
mediator) - is in thermal equilibrium, the virial theorem is not 
satisfied for the mediator, namely T p (T) — T X (T) is non-zero 
for any finite 7. This is a manifestation of quantum entan- 
glement between a particle and a single thermal bath |7J]. The 
difference between the average kinetic and potential energy of 
the mediator can serve as a criterion for the coupling strength 
also in the non-equilibrium case, as shown below. 

Eqs. ( [T2I 1. ( fT6l l. and ( f34t . and the relation p = mx, one can 
obtain 



T = 



and 



hiUQTpS 



t p S = ^2g n [I n {l-fi 2 n )+I nl +I n2 ] (64) 



Tl=l 



T P = JT J2 9n^l [I„(l - £») + I n i + I n %] ■ (65) 



Here An = H n /D, 



xdx 



i[arctan(a) — (7r/2)sign(b)] — ln(|/i n | 2 ) 

2(1 -An) 



(66) 



With a = (/4 r - /4i)/(/-WMra)> & = MnrMni. Mnr = Re(//„) 

and /i„j = Im(/i n ). Finally, 



xdx 



{e x - 1)[X 2 + {ThvjJLn)' 



(67) 



where r h i, 2 = h/k p T lj2 - 

The weak and strong coupling can be defined in terms of 
the effective bath-particle interaction strength 70, where 

7D = D 2 ^{D 2 +loI)- 1 . (68) 

If 7d -C ujq, one can find from dT9b-(L2TT> that 

Mi, 2 ~ =Rw + 7d/2, /i 3 = D - 7 D (69) 

and 

7d 7d 



51,2 ~ T 



2{D 2 



93 



D 2 



(70) 



If 7d *C Wo, the mediator can be described as an oscillator 
with relatively small effective friction. However, this condi- 
tion by itself is not sufficient to guarantee the virial theorem. 
Only if temperatures of both thermal reservoirs are sufficiently 
high, i.e. 

fc B T 1 1 , 2 /^»7 D , (71) 
we get approximately equal steady-state values of T x p : 

ter x w*Br p «i[i7(Ti) + ^(r2)]. (72) 

Here U(T) = hu Q /2 + hui [exp(Huj a / k B T) - is the av- 
erage energy of a quantum oscillator in thermal equilibrium. 
The inequality ( |7TT i can be considered as the usual condition 
of applicability of the Gibbsian statistics, when interaction en- 
ergy between subsystems of a large closed system is small 
with respect to the internal energies of the subsystems l36Tl . 
Thus, as long as the virial theorem is preserved, one can as- 
sign to the mediator a certain temperature T on the basis of 
Eq. (El: U(T) = (l/2)[J7(Ti) + U(T 2 )]. In the high tem- 
perature limit knTi t 2 ^ ft^o, this leads to 



T x w T p w T w (Ti + T 2 )/2, 



(73) 



as is expected. 

For the case of moderate or strong coupling (overdamped 
mediator), T x p acquire 7-dependence and the condition ( l72b 
is not satisfied any more. Fig. |2] shows the monotonic de- 
pendence of the relative energy imbalance AT px = (T p — 
T x ) / (T p +T X ) on the coupling strength. The cut-off parameter 
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FIG. 2: (Color online) Dependence of AT px on the coupling strength 
for different D/uio, Thiuio, and tm^o- (a) D/ujq = 10, Thi^o = 
T"h2^o = 1, (b) D/lu = 1, r h iCJo = T ll2 aj = 1 (c) D/ojo = 10, 
Thiuo = Th2^o = 0.5, and (d) D/uj = 1, T hl uj = 1, t 1i2 ojo = 
0.5. 



can loosely associate with the Debye frequency. The ther- 
mal bath modes with frequencies higher than D are effectively 
decoupled from the mediator and do not play a significant role 
in the processes of thermalization and energy transfer. For this 
reason, when D/cjq decreases, AT px also decreases. This also 
follows from d68b . Indeed, the coupling constant 7 = 7/m is 
renormalized by the factor D 2 / (D 2 + us 2 .) and is effectively 
determined by 70 at a relatively small 70- As one also finds, 
AT px decreases when T\ 2 grow in accordance to d7Tb - d73b . 
Thus, ATp X can be considered as a measure of the coupling 
strength that takes into account all the relevant factors such as 
D 2 /(D 2 +uj 2 ) and values ofTi j2 . 

It is interesting to notice that the 7-dependence of the heat 
conductance d43l ) may possess a maximum. As we found, it 
appears at relatively small ojq and relatively large 7 (strong 
quantum entanglement between the particle and baths) when 
T h uj < 2, h^/k B T > 10, and j/D > 1 as is illustrated 
in Fig. [3] The maximum strength can be characterized by 
the quantity AK = (K max - Koo)/K max , where AT max and 
Koo are the values of K at its maximum and at j/ojq —> 00, 
respectively. As our simulations show, AA' increases when 
Th^o decreases and the value (7/wo)max at which K max is 
achieved shifts toward larger values when D/ljq decreases. It 
is worth to mention that a similar maximum in the dependence 
of the heat current on the system-bath coupling strength was 
found in l24ll for Heisenberg XY spin 1/2 chain. 

Fig. [4] shows the temperature dependences for the normal- 
ized heat conductance at different D and for luot p = 10. The 
shown dependencies are based on the same expression d43l l 
and its high-Z) limit. As our analysis reveals, if D/ujq > 10, 
K{T) essentially coincides with the corresponding D — > 00 
limit d62l ). The straight line region corresponds to the low- 
T limit d5Tb , which is the same for all curves. At large T, 
each curve reaches its classical plateau in accordance to ( 1451 ). 
An unusual feature, which is an additional plateau that ap- 
pears in the intermediate range of r p /rh and D « ojq, is 
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FIG. 3: (Color online) Dependence of the normalized heat conduc- 
tance K/(k-BUio) on the coupling strength at D/ljo = 1 and different 

T h U . 
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FIG. 4: (Color online) Temperature dependence for the normalized 
heat conductance A'/(fcBtJo) at ojot" p = 10 and different values of 
D. Note that r p /r h = kvT/hrf. (a) Dt p = 0.1, (b) Dt p = 1, (c) 
Dr p = 100, and (d) Dr p = 00. 



a generic result valid in a wide range of ljot p . It can be ex- 
plained in the following way. If T is so large that WoTh < 1 
or Tp/rh > ^o T pi we have the "final" classical plateau in 
Fig. 2] described by d47b . On the other hand, if T is so 
small that Dr^ > 1 or T p /Th < Dt p , we have the quan- 
tum regime (straight line in Fig. |4| described by ( BTI ). In 
the intermediate region, when Dt^ < 1 but ojqt^ > 1 or 



Dt p < Tp/Th < ^oT p , one can approximate uj- 
under the integrals in d43T > for n = 1 or 2 and 

rlk B D 2 



K 



87TT D 



Ji ^2 9n+93D 2 J 2 



where 



Jk 



duiuj 



Ml,2 



(74) 



(75) 



/o smh 2 (phu; /2)(D 2 +Lu 2 ) k ' 
Taking into account that g\ + g 2 = —.93, one can rewrite ( |74l i 



8 



o- 



■o- 

K, 



FIG. 5: (Color online) Chain of macroscopic subsystems ("nanopar- 
ticles") interconnected by mediators. Each subsystem indicated by a 
large circle, as well as both thermal reservoirs, correspond to Hamil- 
tonian {3). Other symbols also have the same meaning as in Fig. Q] 
Temperatures Tl and Tr of the left and right thermal reservoirs are 
fixed because their heat capacity is considered infinite. The temper- 
atures T n of the subsystems can slowly vary until the steady-state is 
established. Thermal conductances K n can vary from one connec- 
tion to another. 



8lTT n 



duui 



o smh\f3hio/2)(D 2 +u 2 y 



(76) 



Finally, using that Di\ < 1 and approximating sinh(aj) w x, 
the integral in (|76*| l can be estimated as tt/(Dt£). Due to d70t . 



if 



fc B 5 fc B 7 D 3 j 
-33 = 



8ttt d 



2^ 



(77) 



does not, indeed, depend on temperature forming the plateau 
in Fig. |U The physical origin for the small K and intermedi- 
ate plateau can be explained as follows. If T increases above 
the Debye temperature 6*d = Dh/ks, i.e. Dt^ = Qr>/T < 
1, all the baths modes begin to be excited. At the same time, 
if T is still less than hoj^/k-Q (or t^ujo > 1, which is always 
possible if ljq 2> D), the quantum system cannot be excited 
and, hence, it cannot absorb energy from either bath and trans- 
fer it to the other bath. This will lead to a small value of 
K in (1771) . Moreover, because this situation stays unchanged 
when T changes within the interval 0rj < T < hujQ / k-Q (or 
Dt p < Tp j t h < ujQTp), K must not depend on T notice- 
ably, which is in accordance to d77l ). Also, as one can notice, 
the effective bath-particle interaction strength 7d is small if 
D 3> ujo even if 7 itself is not small. In this case, we have a 
very low decay rate for the /ii 2 modes and, correspondingly, 
very small heat current and heat conductance, again in accor- 
dance to ( T771 i. This result brings our situation close to the one 
with the bound states mentioned at the end of the first part of 
Sec. [TIT1 Indeed, when 7d = 0, we would have non-decaying 
oscillatory contributions to the steady-state, depending on the 
initial conditions, and zero heat current. However, even for a 
case with infinitely small 7d, the unique but infinitely small 
steady-state current will be established after an infinitely large 
time interval. 



VI. FOURIER'S LAW 

The results obtained in previous sections will be applied to 
an extended model that can shed additional light on the long- 
standing problem regarding the origin of Fourier's law lfl4ll . 



The standard definition of a macroscopic body in the state of 
weak non-equilibrium l37ll is that it can be divided into re- 
gions large enough to be considered macroscopic, but small 
enough to be described by a local temperature. In addition, 
these regions or subsystems must be weakly coupled to each 
other. The weak coupling in this context means that the char- 
acteristic time required for the subsystems to come to mutual 
equilibrium is much longer than the time of microscopic re- 
laxation. Here we introduce a model that fills this conceptual 
framework with a microscopic content. 

Fig. [5] illustrates our model. It consists of N macroscopic 
subsystems and two thermal reservoirs (TR) coupled by the 
mediators. Each subsystem and each TR is described by the 
Hamiltonian given by Eq. (O within the framework of the 
Drude-Ullersma model, Eq. (TBI ). Each mediator is described 
by the Hamiltonian and each coupling is described by the 
Hamiltonian (0]i and Eq. (fl4] >. The total system, including the 
TRs, is Hamiltonian with constant total energy. 

Each subsystem and each TR is initially prepared in the 
state of thermal equilibrium and is characterized by a tem- 
perature T n or Tl and Tr, respectively. It means that the re- 
spective correlators have the form of Eq. fl35l l. 

The energy of a given subsystem, which consists of very 
large or infinite number of modes with frequencies Uk = fcA 
with k = 1, 2, 3 ... , is 



U(T) 



(fc B r) 2 

HA 



Here, the divergent zero-point energy term was dropped. Cor- 
respondingly, the heat capacity of a given subsystem 



k B T 

'' hA ■ 



Thus, the thermal reservoirs will be characterized by A — >• 
and, correspondingly, an infinite heat capacity, while the sub- 
systems need to be described by an arbitrarily small, but still 
finite A. The energy difference between two subsystems (as 
in Fig. Q]) with temperatures T\ and T2 is 



AU 



hA 



hA 



On the other hand, the heat current between the two subsys- 
tems can be estimated on the basis of Eqs. d45b and (08]) as 



Jth~— |Ti-T 2 | 



Thus, the characteristic time of mutual equilibration 



"eq 



AU 

Jth 



' hA 



Th 



A" 



(78) 



For the subsystems, we assume that the "Heisenberg" time 
scale 1/A is much larger than any other characteristic times, 
such as r p , r x , or Th 00]. Thus, the system described by the 
Hamiltonian represented by the diagram in Fig. [5] meets the 
conditions described in the first paragraph of this section. 
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One of the main results that we have obtained by solving the 
Hamiltonian (fl} is that the energy flows from higher to lower 
temperature TR. Since all the modes of a given TR are in ther- 
mal equilibrium at the same temperature, the thermodynamic 
relationship between energy and entropy dUi = TdSi is sat- 
isfied for each of them and the total entropy in the steady-state 
increases 



f t =^ 



1 



1 



> 0. 



One can easily show that for the system represented in Fig. [5] 
the total entropy also increases even when the energy currents 
between the subsystems are all different. In the steady-state, 
which corresponds to all energy currents between two neigh- 
boring subsystems being equal, the entropy increases as 



dS 



= |Jth| 



1 1 



> 0. 



These results are valid as long as all modes remain in ther- 
mal equilibrium or close to equilibrium. The coupling con- 
stants given by Eq. ( fT4b are proportional to the infinitesimal 
parameter A, 

C v ,i ~ A 1 / 2 . 

The effect of such coupling on the mediator is finite because 
all modes contribute constructively. The rate of change of the 
correlators ( f35T > is much slower because each of the modes is 
coupled only to the mediator with vanishingly small coupling 
constant. One can see from Eqs. (O - ( fT2b that the rate of 
change of the correlators d35l ) 

d d 

-7^(Pvi(t)Pv>j(t)) ~ —(x m (t)x„,j{t)) ~ A. 

Thus, if we consider the evolution of the system on the time 
scale t, such that 



Tp <C t <§C ^eq; 



(79) 



which is long enough for the microscopic relaxation to take 
place, but short on the macroscopic time scale, all modes will 
remain approximately in thermal equilibrium determined by 
the initial conditions. 

The same argument allows us to extend the solution of the 
Hamiltonian ([TJ shown as a diagram in Fig. Q]to the Hamilto- 
nian that corresponds to Fig. [5] Imagine that after we prepared 
both TRs and all the subsystems in the state of thermal equi- 
librium at the corresponding temperatures we start turning on 
the couplings to mediators one-by-one from left to right. The 
energy flow between two subsystems indexed as n and n + 1 
will be the same as for the Hamiltonian (Q~|i because the effect 
of the subsystem n being already coupled to the subsystem 
n — 1 adds to the energy current a contribution of the order of 
A and therefore negligible on the time scale (1791 . This means 
that the solutions of the Hamiltonian ((TJ can be directly ap- 
plied to the Hamiltonian of the chain shown in Fig. [5] in the 
form of energy conservation condition 



Here the energy currents J n -i,n are given by Eq. d42l and 
this equation is valid for arbitrary values of the initial temper- 
atures. In order to obtain Fourier's law we have to take the 
temperature differences between the neighboring subsystems 
small and express the currents in terms of thermal conduc- 
tances (l43l 

d t U n = K n _i >n {T n ^i — T n ) — K 7hn+ i(T n — T n+ i). (81) 

This equation can be rewritten in the differential form by in- 
troducing a continuous coordinate x, where x ~ nd corre- 
sponds to the locations of the subsystems and d is the distance 
between them. For identical mediators Eq. dSTl i can be recast 

as 

C(T)d t T(x) = K(x - d/2)[T(x - d) - T(x)] - 

K(x + d/2)[T(x)-T(x + d)], (82) 

where C(T) = dU/dT is the heat capacity of a subsystem. 
Then, Eq. (|82l leads to the energy conservation condition with 
Fourier's form of the energy current 



C{T)d t T{x) = d x [ K (x)d x T(x)}, 



(83) 



$t U n J n — 1 ,tj Jn,n 



+ 1- 



(80) 



where C = C/d is the specific heat of the chain, and the 
thermal conductivity k(T) = K(T)d. 

This analysis shows that as long as the energy flow between 
macroscopic subsystems satisfies the condition of entropy in- 
crease, namely, that the energy flows from higher to lower 
temperature, Fourier's law is a straightforward consequence 
of energy conservation. What has been proven in this paper 
is that the dynamics of the Hamiltonian system described by 
the Hamiltonian ([D and its extension shown as a diagram in 
Fig. [5] does indeed lead to the entropy increase. This state- 
ment is predicated on the condition of local thermal equilib- 
rium or near equilibrium for the subsystems that exchange en- 
ergy between themselves. We have stated in our treatment of 
the Hamiltonian (Q3 that this is an initial condition. As such, it 
remains true for the time interval indicated by inequality ( 1791 . 
However, the generally accepted understanding of the slow re- 
laxation processes 0711 in a macroscopic body implies that as 
each subsystem slowly gains or loses energy, it remains close 
to thermal equilibrium with a certain time-dependent tempera- 
ture due to rapid thermalization on the microscopic time scale. 
Thus, the microscopic derivation of Fourier's law in the con- 
text of our model requires not only to prove that in a certain 
limit (smooth temperature variation) the energy current is pro- 
portional to the temperature gradient, which we have done. 
Equal, if not more important, task is to show that the sub- 
systems described by the Drude-Ullersma model are capable 
of self-thermalization when coupled by the harmonic media- 
tor. This proof will require the study of the dynamics of the 
Hamiltonian (Q3 on the time scale given by Eq. ([78). 

The model depicted in Fig. [5]resembles to some extent the 
models analyzed by Michel, Mahler and Gemmer Jit], as well 
as those discussed by Dubi and Di Ventra [,9]. The subsystems 
considered in i||9l] were still microscopic with the heat capac- 
ity of the order of fee- Other publications [10-13] consider the 
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energy transport in chains consisting of spins or harmonic os- 
cillators. In all these cases the goal was to examine the energy 
transport on the "nano-scale". Our approach here is to exam- 
ine a solvable Hamiltonian model of a non-equilibrium system 
that falls within the more traditional, "textbook" framework. 
The subsystems described by the Hamiltonian (O can have 
arbitrarily large heat capacity determined by the infinitesimal 
parameter A. Correspondingly, these subsystems remain in 
thermal equilibrium for the extended period of time, much 
greater than the microscopic relaxation time of the mediator 
and demonstrate that the energy flow leads to entropy increase 
and, hence, to Fourier's law, Eqs. d80l - (l83l ). The mediators 
do not have to be in the state of thermal equilibrium for the 
Eq. ( l83l to be valid. As was discussed in Sec. V, in the 
strong coupling regime the mediator cannot be assigned a cer- 
tain temperature. However, as long as the subsystems remain 
close to thermal equilibrium, Fourier's law is still valid. 

The outstanding question of the quasistatic evolution and 
self-thermalization of such subsystems on the much longer 
time scale, Eq. ( 1781 . will be addressed elsewhere. 



Q '^'^ // q '• [ 

T, T 2 I I T N 

Energy flux 

FIG. 6: (Color online) Diagram of the quantum thermal baths model 
equivalent to the Hamiltonian shown in Fig. [5] Here F n are the 
random forces that appear in Eq. j841 > and the other elements are 
similar to ones from Fig. [5] Notice that the mediators are not coupled 
to each other directly. 



with phonons, their mean free path would be minimum possi- 
ble, equal to "interatomic" distance. 



VII. HEAT CURRENT BETWEEN STM TIP AND 
SUBSTRATE 



A. Quantum thermal baths model 

In the limit of large D, Eq. (fSOb follows from the model 
shown in Fig. [6] The equation of motion for the n's mediator 
is practically the same as Eq. d53t 

mx n + jx n + muj^Xn = F n _i(t) + F n (t), (84) 

where the random force F n describes the effect of the corre- 
sponding subsystem n on the mediator and the correlator of 
the force is determined by the temperature of the respective 
subsystem. Note, that the mediators are not directly coupled 
to each other by "springs", so that Eqs. ([84-b are not a system 
of coupled equations, but simply are n identical equations, 
whose solution for the heat current J n -i n lS given by Eq. 

4SB- 

One can compare this model with the models of self- 
consistent reservoirs ll38l440ll . Since our subsystems are 
macroscopically large, their effect on the mediator can be de- 
scribed in terms of a stochastic force. However, these forces 
are not arbitrarily introduced. We have shown that the exact 
solution of the Hamiltonian dynamics in the limit of large D 
yields the same expression for the energy current as the quan- 
tum thermal baths model. In this sense the quantum thermal 
baths model can be a useful shortcut, but our model is defined 
by the Hamiltonian ([T). It is also important to mention that the 
thermalization of these self-consistent reservoirs is considered 
a given in Refs. I38l - l40ll . In fact, it has to be proven by analy- 
sis of the long term evolution of the Hamiltonian dynamics. 

One should note that in our model the mediators are not 
coupled to each other directly. The energy flows only through 
the macroscopic subsystems. In this sense the chain of media- 
tors operate, by design, in the minimum thermal conductivity 
limit because the neighboring mediators are always uncorre- 
lated. If we were to associate the movement of the mediators 



The results obtained in previous sections can be used to 
clarify an interesting finding reported recently by Altfeder, 
Voevodin, and Roy 12511 . In their experiment, the energy cur- 
rent in high vacuum between the Pt/Ir tip of the scanning tun- 
neling microscope and gold film serving as a substrate was 
found to be anomalously large, exceeding by ten orders of 
magnitude the current given by the blackbody radiation the- 
ory. An interpretation of this phenomenon given in Ref. l25ll 
involves "emission" of phonons by Au surface, facilitated by 
the electric field, and their "tunneling" through the vacuum 
gap to the tip. Here we would like to offer an alternative 
mechanism that seems to be capable to account quantitatively 
for the same effect. 

An important feature of this experiment is that a carbon 
monoxide (CO) molecule was always present between the Pt 
tip and the substrate. From the value of the tunneling conduc- 
tance, the gap between the "last tip atom" and the surface of 
the Au film was inferred to be close to 3 A lT4lTl . The diameter 
a of a CO molecule is close to 3.7 A l42tl . Thus, without the 
CO molecule, the distance between the tip and the substrate 
would be about 7 A. There is also a strong electric field in 
the gap, which is due to the work function difference between 
Pt and Au (w 0.7 eV). The implication of this is that the CO 
molecule must be strongly coupled to both the tip and the Au 
surface. If we consider that this molecule serves as a mediator 
between the two thermal reservoirs (the tip and the substrate, 
similar to the arrangement in Fig. 1), the energy current be- 
tween them can be estimated and compared with the experi- 
mental data. The coupling strength between the mediator and 
the thermal reservoirs determines how strongly damped the 
mediator is. 

There are two vibration modes associated with CO 
molecule attached to the Pt tip l43ll . One of them is the Pt- 
C stretching vibration with frequency 480 cm -1 . The other 
mode is the C-O stretching motion with the frequency ap- 
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proximately 2000 cm -1 . Using the conversion coefficient 1 
K = 0.7 cm -1 , we find that the lowest frequency 480 cm -1 
corresponds to 685 K and the highest frequency 2000 cm -1 
corresponds to 2.8 xlO 3 K. In the experiment, the tip was 
maintained at room temperature, while the temperature of the 
substrate was substantially lower. Thus the higher mode was 
clearly not activated and we should take the frequency of the 
mediator in our model loq rj 480 cm -1 . 

Since the main frequency Swo/^b ~ 685 K is substantially 
greater than the temperature of the bulk of the tip (w 300 K) 
and the temperature of the substrate (in the range 90-210 K), 
the first order approximation of the heat current mediated by 
the CO molecule can be estimated using Eq. ( Bil l 



31,4-2 



Jco 



7t 6 k 



B 7 



(85) 



The expression for the energy current between the tip and the 
substrate due to phonon emission from the "hot" spot on the 
surface was used in Refs. 1251 44|l in the form 



Ja 



m 4 -T 2 4 ). 



(86) 



Here 6*d is the Debye temperature of the substrate. The un- 
derlying physics that leads to these two expressions is very 
different, but the temperature dependence is the same. The 
analysis based on Eq. d86b leads to a good description of the 
experimental results. Thus, the model of CO-mediated heat 
transfer will also give the same results if the prefactors of Eqs. 
( |85| > and (|86l l are equal. This requires that 



7T HujQ 

V2k B 0D 



(87) 



For gold 6*d ~ 165 K, and since hcoo/k-B ~ 685 K we get 



9. 



(88) 



Thus, other things being equal, the model based on CO- 
mediated heat exchange will give the same quantitative results 
as the model of phonon tunneling, provided that the coupling 
between the CO molecule and both the tip and substrate is 
strong enough to make the Pt-C vibrating mode overdamped 
according to Eq. ( f88l >. This is likely the case, considering (as 
mentioned above) that the size of the molecule is comparable 
to the gap between the tip and substrate and also taking into 
account the presence of the strong polarizing electric field. 

The heat current given by Eq. d85l ) is ~ 10 10 times greater 
than that determined by the black-body radiation, Eq. d52l ). 
emitting from the area equal to the cross-section of the CO 
molecule A = 7ra 2 /4, with a ~ 3.7 A, provided that j/ujo ~ 
10. 

Our conclusion is that the anomalously large heat current 
between the STM tip and the substrate can be understood 
just as well as the effect of mediation by the CO molecule. 
There are ways to modify the experiment in order to deter- 
mine which of the two mechanisms is responsible for the ef- 
fect. One is to carry out a similar measurement, but without 



the CO molecule lodged in the tunneling gap. Another op- 
tion is to use the tip and substrate made of the same metal, so 
that there will be no work function potential difference A<I> 
between them. The mechanism of phonon tunneling is based 
on the interaction between the electrically charged tip and its 
electrostatic image. The amount of charge is mainly deter- 
mined by the large electric field E = A$Au/Pt/^> where d 
is the vacuum gap l25ll . In the absence of such large field the 
mechanism of phonon tunneling should be greatly weakened. 



VIII. MINIMUM THERMAL CONDUCTIVITY 

The topic of minimum thermal conductivity can be first 
traced to Einstein's contribution (see Ref. 14511 and references 
therein). In strongly disordered solids the phonon coherence 
length may become of the order of interatomic distance and, 
obviously, cannot be reduced any further. Correspondingly, 
the thermal conductivity reaches its minimum value, at least 
as far as its dependence on such length is concerned. A de- 
tailed treatment of this problem is given in Refs. j46l l47ll . 
The same phenomenon leads to a minimal electric conductiv- 
ity in disordered conductors - the so-called Mott-Ioffe-Regel 
limit 11481 14911 - when the coherence length of the charge car- 
riers become comparable to the interatomic distance. 

Here we would like to add another perspective on this 
matter using the results obtained above and making use of 
the scaling approach previously developed for electron trans- 
port II50h52I1 . It should be emphasized that this section is not 
intended as a comprehensive treatment of this problem, but 
rather as a brief introduction to an alternative approach which 
may be useful in some systems. 

Consider a microscopic block whose edges are along the 
principal axes {x, y, z} of the thermal conductivity tensor. We 
can choose the sides of the block such that, on average, the 
random phase acquired by phonons due to an inelastic inter- 
actions along the way between the two opposite boundaries is 
the same, of the order of 2ir, for all three pairs of the block 
boundaries. This choice for the sides of the block corresponds 
to the definition for the anisotropic phase coherence lengths 
tip,% (i = V: %})> the distances over which phonons lose 
phase coherence. 

Let K v: i be the thermal conductance of such a phase co- 
herent volume (PCV), so that the energy current through this 
block 



J<p,' 



K Vti 5Ti. 



(89) 



Here STi is the temperature difference between the opposite 
edges of the PCV block. Notice that the notion of a tempera- 
ture difference cannot be introduced for distances smaller than 
the phase coherence length. In order to express the macro- 
scopic anisotropic thermal conductivity m in terms of the con- 
ductance of the PCV, consider a macroscopic block with sizes 
{L x , L y ,L z } obtained by fitting together N 3 phase coherent 
volumes, so that L X /£ V)X = L y /£ v , y = L z /£ VtZ = N > 1. 
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By virtue of Fourier's law, the heat current in the 
x— direction through this macroscopic volume is 



J x = k, A/'.— r— 



(90) 



Here we assume a linear temperature variation across the 
block. On the other hand, the heat currents through the in- 
dividual PCVs combine linearly, so that 

J x =K ip , x 5T x N 2 . (91) 

Taking into account that 5T X = AT X /N, we obtain 



K v>x - v f . (92) 



there are no propagating phonons. Instead, every atom is acted 
upon by the non-equilibrium environment and the mechanism 
of this interaction can be described by the Hamiltonian given 
by Eq. (see also Fig. 1). Thus, the conductance K v of the 
PCV containing one atom can be well described by our model 
in which a single oscillator mediates the energy exchange be- 
tween two thermal reservoirs. Then, by virtue of Eq. d92l the 
thermal conductivity is given by 



(94) 



<p< x p p 

To put this result into a different perspective, the PCV is the 
minimal size block for which one can introduce the notion of 
thermal conductivity. 

One can argue that the thermal conductance of the PCV 
of phonons in an anisotropic medium is isotropic, similar to 
the electric conductance in anisotropic metals 15311 . namely, 
K v . x — K ViV = -K^z. If this assertion were true, the 
anisotropy is defined by the following relationship: 



where £q is a constant of the order of the interatomic distance 
and Ku, = K is given by the general expression, Eq. j43l , 
and its limiting cases such as Eq. d48b . 

Let us consider the classical limit of high temperatures, Eq. 
(l48l i. when Eq. d94l l takes the form 



4t, 



&B7 

44 ■ 



(95) 



This has to be compared with another expression for the MTC 
which is based on the atomic density n and elastic constants 



ft a 



£ 2 

~p~2 ' 

<p,v 



K>3 



e 2 



K min = 0Ak B n 2/3 (vi + 2v t ). 



(96) 



(93) 



This relationship is in agreement with the quasiclassical re- 
sults obtained from the kinetic equations /tj ~ cUjAj, where 
c is the specific heat, v,i is the anisotropic speed of sound, 
and Xi is the phonon mean free path. Considering that the 
mean free path is similar to the phase coherence length, so 



Here v\ and u t are the longitudinal and transverse speeds of 
sound, respectively. This is the sum of the contributions of 
the three spatial degrees of freedom. Taking into account that 



,2/3 



1 



1 



that A, 



ViT v , we get Ki 



the relaxation (decoherence) time t v and the specific heat are 
scalars, the relationship (l93l follows. 

Eq. ( |93l is most useful when applied to a system where at 
least one of the coherence lengths is temperature independent 
constant. It may be a highly disordered crystal in which the 
decoherence takes place over interatomic distances, or a lay- 
ered structure in which the coherence length in one direction 
is fixed by the size of the layer. Further discussion of the con- 
sequences of Eq. ( |93l would be far outside the scope of this 
paper. 

Now we will return to the case of highly disordered sub- 
stances in which the coherence lengths in all directions are 
of the order of interatomic distances and do not change with 
temperature. There are numerous example of such substances 
where the minimum thermal conductivity is reached at tem- 
peratures above 30 K j45ll . In vitreous, silica- and germania- 
based glasses, the mean free path (or the phase coherence 
length) approaches the interatomic distance at T > 100 K. 
For example, the phonon mean free path for amorphous sele- 
nium at T > 50 K is temperature independent and equal to 
5 x 10~ 8 cm, which corresponds app roximately to the inter- 
atomic distance in this substance 15411 . 

In the minimum thermal conductivity (MTC) regime the 
movements of the neighboring atoms are incoherent, so that 



c11Jt v . Since both K rai „ can be estimated as 



£o to ' 



(97) 



Hereafter we will drop the numerical prefactors. The charac- 
teristic time scale £q /v av is the time of flight of a phonon over 
the interatomic distance. The incoherence of the neighboring 
atoms means that the phonons lose their coherence over this 
time interval. This is exactly the meaning of the decoherence 



time t v - £ /v 



In our model the relaxation time r p is de- 



termined by the strength of the coupling between the oscillator 
and the thermal baths. It is the relaxation time of the momen- 
tum of the oscillator. Thus, in the MTC regime t v and t p are 
equivalent quantities, 



(98) 



and the expressions for the MTC given by Eqs. d95l ) and d96l ) 
are qualitatively and even quantitatively similar. 

Moreover, since the standard definition of the Debye fre- 
quency is wd = k^O-D/h ~ v av /£o, we see that in the MTC 
regime the relaxation time in our model must be 



7 — w D - 



(99) 
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The frequency loq is the highest frequency associated with 
atomic vibration and it must be also of the order of Debye fre- 
quency, ujq ~ v av /£o ~ w d- The last parameter of the model 
is the cut-off frequency D, which defines the maximum fre- 
quency of the modes of the thermal baths that are coupled to 
the mediator. In the context of solid substances this cut-off 
also must be of the order of the Debye frequency. Thus, we 
come to conclusion that the range of parameters within which 
the Hamiltonian model given by Eq. (1) is applicable to the 
description of the minimum conductivity regime is rather nar- 
row and is given by 

uj ~ 7 ~ D ~ uj d . (100) 

As an example, we can take the data from Ref. ifsill for 
selenium at temperatures above 100 K. The value of thermal 
conductivity K m i n ~ 0.5x 10~ 2 WK _1 cm _1 . Comparing this 
value with Eq. (l95l l, and taking into account that the character- 
istic interatomic distance £q w 5 x 10~ 8 cm, we find 7 ~ 10 13 
s _1 . The Debye temperature for selenium is 6*d ~ 250 K, so 
that wd ~ 3 x 10 13 s _1 and the condition given by Eq. ( |99] l 
is satisfied. Thus, in a highly disordered substance the oscilla- 
tor enclosed inside the PCV is rather overdamped, ujqt p > 1. 
Although this is only an estimate, we use it in order to illus- 
trate the potential applications of our model. A more detailed 
comparison with the experimental data needs to involve the 
specific heat also calculated within the framework of the same 
model. 



IX. JOSEPHSON JUNCTIONS 

Finally, we can mention that for some potential applications 
of our model, such as the Josephson junctions, all model's 
parameters are already experimentally known. This enables 
making valuable predictions about the physical behavior of 
the corresponding system. In the case of the Josephson junc- 
tion, the particle's coordinate x in the Langevin equation is 
substituted by <fi which is the phase difference between the 
wave functions describing the state of condensate of Cooper 
pairs in the contacting superconductors kept at different tem- 
peratures. Here luq is the plasma frequency and 7 = 7/m = 
1/RC with R and C are the junction resistance and capaci- 
tance, respectively |[55ll . Characteristic values for ujq can vary 
between 10 10 s _1 and 10 14 s~\ depending on the current den- 
sity. For tunnel junctions, when 7/wo is usually in the range 
0.001 - 0.1, the damping is weak. For junctions with non- 
tunneling conductivity and in the form of point contacts or 
thin-film microbridges l56ll . 7/t^o >> 1 and the damping is 
large. 

X. SUMMARY 

In conclusion, we have considered the heat transport be- 
tween two thermal reservoirs mediated by a quantum parti- 
cle using the generalized quantum Langevin equation. Both 



thermal reservoirs are described as ensembles of harmonic 
modes using the Drude-Ullersma model and the mediator is 
also treated in the harmonic approximation. The expressions 
obtained for the heat current and thermal conductance are 
valid for arbitrary coupling strength between the mediator and 
the reservoirs. The cutoff frequency, which characterizes the 
thermal reservoirs, can also be arbitrary. The obtained results 
are analyzed for different temperatures regimes and different 
strengths of the coupling parameter. The dependence of the 
thermal conductance on the coupling strength shows a maxi- 
mum and the temperature dependence of this quantity reveals 
a plateau at intermediate temperatures, similar to the classical 
plateau that corresponds to the high-temperature limit. 

The results are applied to a model of a chain made out of 
macroscopically large, but finite subsystems, each described 
by the Drude-Ullersma model. These subsystems are cou- 
pled to each other through a quantum mediator. As long as 
the subsystems are large enough, so that their energy changes 
slowly in comparison with the relaxation rate of the media- 
tor's energy, Fourier's law follows as a differential form of 
the energy continuity equation. It is important to notice that 
at no point this derivation relies on any assumptions outside 
the framework of the Drude-Ullersma model. Thus, it may be 
considered as one of the few examples of rigorous derivation 
of Fourier's law from the first principles, at least on the time 
scale that leaves the modes of the thermal baths in thermal 
equilibrium. 

We have applied our results to explain the observed anoma- 
lously large heat flux between STM tip and substrate. Our 
conclusion is that the effect is due to the mediating role of 
the CO molecule placed in the tip-substrate gap. We also out- 
lined the approach by which our model can be applied in or- 
der to understand thermal conductivity of highly disordered 
substances - the minimum thermal conductivity - and to the 
non-equilibrium Josephson junction. 
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APPENDIX 

In the second way, at t < 0, the dynamical variables x-n {t) 
and P2i{t) of the second bath are determined by relations to 
(l32l and (l33l l. while xu(t) and xu(t), which now incorporate 
the quantum system, are determined as 



xu(t) = J2 \ L H . ^«^ lkt + (A101) 
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and Pii(t) = muX\i{t). Here ef are orthonormal eigenvec- 
tors Il7J]: 



(a) ^-(b) _ ~ T-(a) 



(A109) 



k I D 2 + v\ 2Aoji sin 0(w fc ) 

Taking into account ( f34b and ( lAlOU , one finds 

ft 

2 x /m 1?; m lj ^ z/ fe 



(A102) 



(^(0)^(0)) 



J coth-i— -^,(A103) 



Using i57|, the above summations can be carried out accu- 
rately and the result is 



1 w 2 D 2 + v 2 sm<t>{v) 



(pii(o)pij(o)> = — 2^ " ke * e * coth 



2 



,(A104) 



and (pii(0)xij(0) + xij(0)pii(0)) = as before. Using these 
relations, the first two sums in d28| i can be recast into 



— 2^ COs(Wit)ji = 3 2^ AUk D + ^ 



x sin 2 ^ coth ^ F[ a \v k ) 4 a \v k ) (A105) 



and 



1 v - ^ n ■ , , v (6) 4ft 7l D 2 ^ £> 
2m f— ^ 7r J m ^ i/fc 



x sin 2 <fk coth V fc )jf (A106) 



In dAT05b and dAlQoI 



Aw 2 cos 



(a) 



n=l,2,3; 



^(ija + w a)( w a-i^) 



g n Au>i [/j, n sin corf — uii cos Wji] 

=1 ( / 4 + u,f)(£2 +w 2 )(w 2_ z , 2) 



(A107) 



,(A108) 



^ a) (-) 



3 



2 L> 2 



i/ 2 ' /i 2 

n=l r 



sin(Vi + </>(V)) cos(ist + (f>(v)) 



sin0(i^) 



sinffl f 



(Alll) 



and J 7 } ' (i/), J 7 ^ (V) are determined from ( IA109I ). In the de- 
rived expressions we disregarded all contributions that are ex- 
ponentially decaying in time. As in the first way, after substi- 
tuting mt0 < IA105b and (IA106b contributions that con- 
tain the product sin(i4 + 4>{v)) cos(vt + 4>(v)) cancel each 
other, and the other time-dependent terms will be proportional 
to sin 2 (Vi + 0(V)) + cos 2 + = 1. As is also clear, the 
coefficient in the product T^'^ T^'^ is inverse proportional 
to sin 2 4>k and is canceled by similar factors in ( 1A105I ) and 
dA106l >. This eliminates the dependence on initial conditions 
related to whether the central particle was initially connected 
or not to the first bath. These observations prove explicitly 
the existence of the steady-state in the presented mode and its 
uniqueness. Finally, replacing the summation over k by the 
integral results in the same expression d40b for (V l ,} < -- ab \ 
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